{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Cell Segmentation</h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Segmentation of cells in fluorescent microscopy is a relatively common image\n",
    "characterisation task with variations that are dependent on the specifics of fluorescent\n",
    "markers for a given experiment. A typical procedure might include\n",
    "\n",
    "1. Histogram-based threshold estimation to produce a binary image.\n",
    "1. Cell splitting (separating touching cells) using distance transforms and a watershed transform.\n",
    "1. Refinement of initial segmentation using information from other channels.\n",
    "1. Cell counting/characterisation.\n",
    "\n",
    "This example demonstrates the procedure on a 3 channel fluorescent microscopy image. The blue channel\n",
    "is a DNA marker (DAPI) that stains all cells, the red channel is a marker of cell death (Ph3)\n",
    "while the green channel is a marker of cell proliferation (Ki67). A typical experiment might count the\n",
    "number of cells and measure size in the different states, where states are determined by presence\n",
    "of Ph3 and Ki67, various times after treatment with a drug candidate.\n",
    "\n",
    "## Acknowledgements\n",
    "\n",
    "The image used in this example is part of the data distributed with the [Fiji training notes](http://imagej.net/User_Guides) by C. Nowell and was contributed by Steve Williams, Peter MacCallum Cancer Centre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cell segmentation and splitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram-based threshold estimation is performed by the <tt>segChannel</tt> function, listed below.\n",
    "It applies light smoothing followed by the Li\n",
    "threshold estimator, one of a range of threshold estimation options available\n",
    "in SimpleITK. A cell splitting procedure using\n",
    "distance transforms and a marker-based watershed (implemented by <tt>segBlobs</tt>, also listed below) was then applied to\n",
    "the resulting mask. Distance transforms replace each foreground pixel with the distance to the\n",
    "closest background pixel, producing a cone-shaped brightness profile for each circular object. Touching\n",
    "cells can then be separated using the peaks of the cones as markers in a watershed transform.\n",
    "A marker image is created by identifying peaks in the distance transform and applying a connected-component labelling.\n",
    "\n",
    "The inverted distance transform is used as the control image for the watershed transform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and display\n",
    "\n",
    "Microscopes use many variants of the tiff format. This one is recognised as 3D by the SimpleITK readers so we extract\n",
    "slices and recompose as a color image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(SimpleITK)\n",
    "## set up viewing tools\n",
    "source(\"viewing.R\")\n",
    "\n",
    "# Utility method that either downloads data from the Girder repository or\n",
    "# if already downloaded returns the file name for reading from disk (cached data).\n",
    "source(\"downloaddata.R\")\n",
    "\n",
    "# this is to do with size of display in Jupyter notebooks\n",
    "if (!exists(\"default.options\")) \n",
    "{\n",
    "default.options <- options()\n",
    "}\n",
    "\n",
    "cntrl <- ReadImage(fetch_data(\"Control.tif\"))\n",
    "## Extract the channels\n",
    "red <- cntrl[ , , 1]\n",
    "green <- cntrl[ , , 2]\n",
    "blue <- cntrl[ , , 3]\n",
    "\n",
    "cntrl.colour <- Compose(red, green, blue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "show_inline(cntrl.colour, Dwidth=grid::unit(15, \"cm\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the functions that do segmentation and blob splitting for a channel (i.e. separately for red,green blue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "segChannel <- function(dapi, dtsmooth=3, osmooth=0.5)\n",
    "{\n",
    "  # Smoothing\n",
    "  dapi.smooth <- SmoothingRecursiveGaussian(dapi, osmooth)\n",
    "  # A thresholding filter - note the class/method interface\n",
    "  th <- LiThresholdImageFilter()\n",
    "  th$SetOutsideValue(1)\n",
    "  th$SetInsideValue(0)\n",
    "  B <- th$Execute(dapi.smooth)\n",
    "  # Call blob splitting with the thresholded image\n",
    "  g <- splitBlobs(B, dtsmooth)\n",
    "  return(list(thresh=B, labels=g$labels, peaks=g$peaks, dist=g$dist))\n",
    "}\n",
    "\n",
    "splitBlobs <- function(mask, smooth=1)\n",
    "{\n",
    "  # Distance transform - replaces each voxel\n",
    "  # in a binary image with the distance to the nearest\n",
    "  # voxel of the other class. Circular objects\n",
    "  # end up with a conical brightness profile, with\n",
    "  # the brightest point in the center.\n",
    "  DT <- DanielssonDistanceMapImageFilter()\n",
    "  DT$UseImageSpacingOn()\n",
    "  distim <- DT$Execute(!mask)\n",
    "  # Smooth the distance transform to avoid peaks being\n",
    "  # broken into pieces.\n",
    "  distimS <- SmoothingRecursiveGaussian(distim, smooth, TRUE)\n",
    "  distim <- distimS * Cast(distim > 0, 'sitkFloat32')\n",
    "  # Find the peaks of the distance transform.\n",
    "  peakF <- RegionalMaximaImageFilter()\n",
    "  peakF$SetForegroundValue(1)\n",
    "  peakF$FullyConnectedOn()\n",
    "  peaks <- peakF$Execute(distim)\n",
    "  # Label the peaks to use as markers in the watershed transform.\n",
    "  markers <- ConnectedComponent(peaks, TRUE)\n",
    "  # Apply the watershed transform from markers to the inverted distance\n",
    "  # transform.\n",
    "  WS <- MorphologicalWatershedFromMarkers(-1 * distim, markers, TRUE, TRUE)\n",
    "  # Mask the result of watershed (which labels every pixel) with the nonzero\n",
    "  # parts of the distance transform.\n",
    "  WS <- WS * Cast(distim > 0, WS$GetPixelID())\n",
    "  return(list(labels=WS, dist=distimS, peaks=peaks))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Segment each channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dapi.cells <- segChannel(blue, 3)\n",
    "ph3.cells <- segChannel(red, 3)\n",
    "Ki67.cells <- segChannel(green, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DAPI channel provides consistent staining, while the other stains may only occupy parts of a nucleus. We therefore combine DAPI information with Ph3 and Ki67 to produce good segmentations of cells with those markers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a mask of DAPI stain - cells are likely to be reliably segmented.\n",
    "dapi.mask <- dapi.cells$labels !=0\n",
    "# Mask of cells from other channels, which are likely to be less reliable.\n",
    "ph3.markers <- ph3.cells$thresh * dapi.mask\n",
    "Ki67.markers <- Ki67.cells$thresh * dapi.mask\n",
    "# Perform a geodesic reconstruction using the unreliable channels as seeds.\n",
    "ph3.recon <- BinaryReconstructionByDilation(ph3.markers, dapi.mask)\n",
    "Ki67.recon <- BinaryReconstructionByDilation(Ki67.markers, dapi.mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we view the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sx <- 1:550\n",
    "sy <- 1450:2000\n",
    "r1 <- red[sx, sy]\n",
    "g1 <- green[sx, sy]\n",
    "b1 <- blue[sx, sy]\n",
    "colsub <- Compose(r1, g1, b1)\n",
    "dapisub <- dapi.cells$thresh[sx, sy] == 0\n",
    "dapisplit <- dapi.cells$labels[sx, sy] == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A subset of the original - note speckled pattern of red stain in some cells\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "show_inline(colsub, pad=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Segmentation of DAPI channel without splitting - note touching cells on mid right that get separated by splitting process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "show_inline(dapisub, pad=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "show_inline(dapisplit, pad=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets check the segmentation of the Ph3 (red) channel. Note that the simple segmentation does not always include complete cells (see lower right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "show_inline(ph3.cells$thresh[sx, sy]==0, pad=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After geodesic reconstruction the incomplete cells match the DAPI channel segmentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ph3sub <- ph3.recon[sx, sy]==0\n",
    "show_inline(ph3sub, pad=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Characterization and counting\n",
    "\n",
    "Image segmentations can lead to quantitative measures such as counts and shape statistics\n",
    "(e.g., area, perimeter etc). Such measures can be biased by edge effects, so it is useful to\n",
    "know whether the objects are touching the image edge. The classes used for these steps in\n",
    "SimpleITK are <tt>ConnectedComponentImageFilter</tt> and <tt>LabelShapeStatisticsImageFilter</tt>.\n",
    "The former produces a _labelled_ image, in which each binary connected component is given\n",
    "a different integer voxel value. Label images are used in many segmentation contexts, including\n",
    "the cell splitting function illustrated earlier. The latter produces shape measures per\n",
    "connected component. The function below illustrates extraction of centroids, areas and\n",
    "edge touching measures.\n",
    "\n",
    "Cell counts are also available from the table dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function to extract the relevant statistics from the labelled images\n",
    "getCellStats <- function(labelledim)\n",
    "{\n",
    "  # create a statistics filter to measure characteristics of\n",
    "  # each labelled object\n",
    "  StatsFilt <- LabelShapeStatisticsImageFilter()\n",
    "  StatsFilt$Execute(labelledim)\n",
    "\n",
    "  objs <- StatsFilt$GetNumberOfLabels()\n",
    "  ## create vectors of each measure\n",
    "  areas <- sapply(1:objs, StatsFilt$GetPhysicalSize)\n",
    "  boundarycontact <- sapply(1:objs, StatsFilt$GetNumberOfPixelsOnBorder)\n",
    "  centroid <- t(sapply(1:objs, StatsFilt$GetCentroid))\n",
    "  # place results in a data frame\n",
    "  result <- data.frame(Area=areas, TouchingImageBoundary=boundarycontact,\n",
    "                       Cx=centroid[, 1], Cy=centroid[, 2])\n",
    "  return(result)\n",
    "}\n",
    "## Label the cell masks\n",
    "ph3.recon.labelled <- ConnectedComponent(ph3.recon)\n",
    "Ki67.recon.labelled <- ConnectedComponent(Ki67.recon)\n",
    "## Collect the measures\n",
    "dapistats <- getCellStats(dapi.cells$labels)\n",
    "ph3stats <- getCellStats(ph3.recon.labelled)\n",
    "ki67stats <- getCellStats(Ki67.recon.labelled)\n",
    "## begin creating a data frame for plotting\n",
    "dapistats$Stain <- \"dapi\"\n",
    "ph3stats$Stain <- \"ph3\"\n",
    "ki67stats$Stain <- \"ki67\"\n",
    "# Put the data frames together\n",
    "cellstats <- rbind(dapistats, ph3stats, ki67stats)\n",
    "cellstats$Stain <- factor(cellstats$Stain)\n",
    "# Remove cells touching the image boundary\n",
    "cellstats.no.boundary <- subset(cellstats, TouchingImageBoundary == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the data has been collected it can be used for plotting, statistical tests, etc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Reset the plot options after dealing with images.\n",
    "options(default.options)\n",
    "library(ggplot2)\n",
    "ggplot(cellstats.no.boundary, aes(x=Area, group=Stain, colour=Stain, fill=Stain)) +\n",
    "    geom_histogram(position=\"identity\", alpha=0.4, bins=30) + ylab(\"Cell count\") + xlab(\"Nucleus area\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
